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Abstract — Max-plus based methods have been recently de- 
veloped to approximate the value function of possibly high 
dimensional optimal control problems. A critical step of these 
methods consists in approximating a function by a supremum 
of a small number of functions (max-plus "basis functions") 
taken from a prescribed dictionary. We study several variants 
of this approximation problem, which we show to be continuous 
versions of the facility location and A:-center combinatorial 
optimization problems, in which the connection costs arise 
from a Bregman distance. We give theoretical error estimates, 
quantifying the number of basis functions needed to reach a 
prescribed accuracy. We derive from our approach a refinement 
of the curse of dimensionality free method introduced previ- 
ously by McEneaney, with a higher accuracy for a comparable 
computational cost. 

I. Introduction 

Dynamic programming is one of the main approaches 
to optimal control. It leads to solving Hamihon-Jacobi- 
Bellman (HJB) partial differential equations. Several tech- 
niques have been proposed in the literature to solve this 
problem. We mention, for example, finite difference schemes 
and the method of the vanishing viscosity [CL84], the 
antidiffusive schemes for advection [BZ07], the so-called 
discrete dynamic programming method or semi-Lagrangian 
method [CD83], [Fal87], [FF94], [CFF04]. Unhke alternative 
approaches based on the maximum principles or on direct 
methods, dynamic programming based methods are guaran- 
teed to give the global optimum of the problem. However, 
they suffer from the curse of dimensionality, meaning that 
the execution time grows exponentially with the dimension 
of the state space. 

Recently, a new class of methods has been developed after 
the work of Fleming and McEneaney [FMOO], see in partic- 
ular [McE07], [AGL08], [MDG08]. These methods all rely 
on max-plus algebra. Their common idea is to approximate 
the value function by a supremum of finitely many "basis 
functions". They exploit the max-plus linearity of the Lax- 
Oleinik semi-group (evolution semi-group of the HJB partial 
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differential equations associated to a deterministic optimal 
control problem). One of these methods, developed by Akian 
and al. [AGL08], was shown to be a max-plus analogue of 
the (Petrov-Galerkin) finite element method. In particular, 
the global error of the method can be estimated in terms of 
certain elementary projection errors, as in the case of usual 
finite elements. 

Among the max-plus methods, the curse of dimension- 
ality reduction of McEneaney [McE07] (see also [MKIO], 
[MDG08], [SSMIO]) appears to be of special interest. In 
its original form, it applies to an optimal switching prob- 
lem involving m linear quadratic models: it approximates 
the solution by a supremum of quadratic functions which 
are obtained by solving Riccati equations. The theoretical 
analysis of the method [MKIO] shows that the growth of 
the execution time is only polynomial as the dimension 
grows, keeping all the other parameters fixed. However, the 
bound of [MKIO] still grows exponentially as the required 
accuracy tends to zero, hence the curse of dimensionality 
is replaced by a "curse of complexity" [MDG08]. However, 
the complexity of the method can be considerably reduced 
in practice by incorporating a pruning algorithm, which 
eliminates on the fly the redundant basis functions produced 
by the algorithm. In this way, high dimensional instances 
(with state dimensions from 6 to 15) inaccessible by other 
methods could be solved [MDG08], [SSMIO]. 

This raises the question to understand why, and to what 
extent, max-plus techniques can attenuate the curse of di- 
mensionality: this is the object of the present paper 

After a brief review of max-plus based methods (Sec- 
tion [ll|i we establish in Section III a negative result, concern- 
ing the family of max-plus methods based on c-semiconvex 
transforms, developed by Fleming and McEneaney [FMOO] 
and Akian et al. [AGL08]. In these methods, the function is 
approximated by a supremum of quadratic forms all of which 
have the same hessian. Then, Theorem |3 . 3 1 below shows that 
the number of max-plus basis functions necessary to reach 
an accuracy of e is at least of order e^'^/^, meaning that the 
curse of dimensionality is inherent to all of these methods 
(the order e^''/^ is optimal, it is reached in particular by the 
max-plus finite elements of [AGL08], with P2 finite elements 
and Pi or P2 test functions [Lak07]). The proof, which 



is sketched in Section IV relies on results concerning the 
approximation of smooth convex bodies [Gru93], [Gru07]. 

However, this theoretical negative result is contrasted by 
the experimental efficiency of pruning in the dimensionality 
free method of [McE07], which often gives approxima- 
tions of an acceptable accuracy for a modest amount of 
basis functions. Therefore, we focus our attention on the 
algorithmic aspects of the pruning problem in the rest of 
the paper In Section [V] we present a primal variant of 
the method, which avoids the use of dual representations: 
in the absence of pruning, it is equivalent to the original 
method, but we shall see that it leads to a more efficient 



pruning. Next, we show in Section VI that the optimal 
pruning problem can be formulated as a continuous version 
of the A:-median or A:-center problem, depending on the choice 
of the norm. The discrete versions of these problems are 
NP-hard. Hence, we propose several heuristics (combining 
facility location heuristics and Shor SDP relaxation scheme). 
Experimental results are given in Section |VII| They show 
that by combining the primal version of the method with 
improved pruning algorithms, a higher accuracy is reached 
for a similar running time, by comparison with [McE07], 
[MDG08]. 

II. Max-plus numerical methods to solve 

OPTIMAL CONTROL PROBLEMS 
A. The Lax-Oleinik semi-group 

We consider the optimal control problem 

v(x,r) :=sup r£(x(5),u(i))t/i + 0(x(r)) ; (1) 
Jo 

x(s) =/(x(5),u(i)), x(0)=x, x(5) ex,u(.?) e?/ . (2) 

Here, X C W' is the set of states, U C M'" is the set of 
actions, T denotes the horizon, the initial condition x G X, the 
Lagrangian £ : X x U — )• M, the terminal reward (/) : M — )• M, 
and the dynamics f -.X xU are given. The supremum 

is taken over all the control functions u and system trajec- 
tories X satisfying (|2]i, and v is the value function. We will 
assume here for simplicity that the set X is invariant by the 
dynamics Q for all choices of the control function u. Under 
certain regularity assumptions, it is known that v{x,t) is the 
solution of the Hamilton-Jacobi equation: 

-^+//(x,^) =0, y{x,t)eXx{0,T], (3) 
ot ax 

with initial condition: 



v(x,O) = 0(x), VxeX 



(4) 



Let (57')7'>o be the Lax-Oleinik semi-group, i.e., the evolu- 
tion semi-group of the Hamilton-Jacobi equation. For every 
horizon T, St is a map which associates to the terminal 
reward the value function := v{x,T) on horizon 

T. By semi-group, we mean that St+s — St o Ss for all 
t,s > 0. Recall that the max-plus semiring, lR„,ai, is the set 
MU {— oo}, equipped with the addition {a,b) i— max{a,b) and 
the multiplication {a,b) i—?' a + b. For all maps f,g from X 
to Rmitx and X e K„,ax, we denote by / Vg the map such that 



{f\/g){x) — max{f{x),g{x)) and by X +f the map such that 
(A +f){x) = A + f{x). It is known that the semi-group St is 
max-plus linear, i.e.. 



Sr [/ V g] = St [/] V St [g] , S,[X+g] = X+St [g] 



(5) 



We shall see that the max-plus basis method exploit these 
properties to solve the optimal control problem Q. 

B. Max-plus linear spaces 

A set W of functions Kmax is a max-plus linear 

space if for all 0i , 02 € >^ and A G M, the functions 0i V 02 
and A + 01 belong to A max-plus linear space W is 
(conditionally) complete if the pointwise supremum of any 
family of functions of W that is bounded from above by an 
element of ^ is finite. 

Let ^ he a set of functions M'' — > M (max-plus basis 
functions). The complete max-plus (linear) space span^ of 
functions generated by ^ is defined to be the set of arbitrary 
linear combinations of elements of in the max-plus sense, 
so that an element of span^ reads sup^.^ag{a{w) + w) 
for some family {a{w))„^,^ of elements of Mmax- The (non 
complete) space span^ is defined in a similar way, but the 
linear combination must now involve a finite family, meaning 
that a{w) = — oo for all but finitely many values of w G We 
refer the reader to [LMSOl], [CGQ04], [McE06] for more 
background on max-plus linear spaces. 

Several choices of basis functions have been considered 
in the literature. Following [FMOO] and [AGL08], we will 
consider here a set 3§ consisting of the basis functions of 
the form 

c , 



Wp{x) 



\x\ +p X, 



I. 'I 



(6) 



where c is a fixed real constant. Hence, an element of span^ 
can be written as 



0(x) = sup [--\x 



\''+p'^x + a{p)], \/xe'. 



(7) 



for some function a 



Recall that a function 



is c-semiconvex if and only if the function x i— > ^{x) + ^\x\ 
is convex. Then, it follows from Legendre-Fenchel duality 
that the space span^ coincides with the space of c- 
semiconvex (lower semicontinuous) functions. 

When c 7^ 0, it is sometimes more convenient to write the 
expansion (|7| in the form 

0(x) sup[-^|z-x|2+fl'(z)], VxGM'' . 

One can pass from one representation to the other by the 
change of variable p ~ cz. 

If W is a complete max-plus linear space of functions 
R'' — > Kmax, and if is any function W' — > Mmax, the 
projection of onto W is defined to be 



^#-(0) := max{v/ G #' | < 0} 



(8) 



(by writing max, we mean that the supremum element of the 
set under consideration belongs to this set, which follows 
from the completeness of 1^). 



All the previous definitions can be dualized, replacing max 
by min, and — oo by In particular, a complete min-plus 
linear space is a set ^ of functions M'' M.Ll{+°°} such 
that — iF := {~w \ w G 3f} is a complete max-plus linear 
space. Then, we define the dual projector by 

P'^'i^) := min{v/ e | > 0} , 

for all functions (j) -.W' ^ MU {+°o}. 

C. Max-plus basis methods 

All these methods approximate the value function at time t 
by a finite max-plus linear combination of max-plus basis 
functions, i.e., 

v(x,0^vUx) = Vt,(A?+wKx)), 

where V/,w'(x) e Typically, w]{x) — ~f + pj^ 
(see [FMOO]), so that the previous approximate represen- 
tation is nothing but a discretization of the semiconvex 
representation (|7]i of v{x,t). 

Then, the coefficients A/ and the functions w\{x) need to 
be inductively determined. Let us fix a time discretization 
step T, such that T = NT for some integer A^. Using the 
semi-group property, we get 

v(-,f + T)=5,[v(-,0], f-0,l,...,r-T . (9) 

We require the max-plus basis approximation v'/^ of v(-,f) to 
satisfy the analogous relation, at least approximately: 

v;+^^5.[v4] = vti(A/+5.K]), t = o,...,T~T (10) 

The various methods differ in the way they address the two 
subproblems, 

1) for all w', replace 5t[w'] by an easily computable 
(accurate enough) approximation 5t(wJ); 

2) Project each 5t(w^), or v|Lj(A/ to the space 
of finite max-plus linear combinations of basis func- 
tions. 

The first subproblem is the simplest one: computing 5t[w5] 
is equivalent to solving an optimal control problem, but the 
horizon T is small, and the terminal reward w'j (typically a 
quadratic function) has a regularizing and a "concavifying" 
effect, which implies that the global optimum can be accu- 
rately approached (by reduction to a convex programming 
problem), leading to various approximations with a consis- 
tency error of 0{t''), with r — 3/2,2, or sometimes better, 
depending on the scheme, see [McE06], [AGL08], [Lak07]. 

The second step is the critical one, in particular, the 
accuracy of the method is limited by the projection er- 
ror [AGL08], i.e., the maximal distance between every 
function ^^[vy or 5't[vv5] and its best approximation by a 
max-plus linear combination of basis functions. In a number 
of methods, including the original one [FMOO], the approx- 
imation v'l^ is such that, assuming that the approximation 
5t(w;) of 5t(w-) is exact, 

v;,<v(-,f) (11) 



so that v'l^ < P-//rv{-,t) where Wt = spanjwj | 1 < / < qt}, and 
P-gr-^ is defined as in ([S). Then, the approximation error in the 
Lp norm, gp := |jv(-,f ) ~ ^/Jlp satisfies 

ep>\\v{-,t)-Pyy,v{-,t)\\p (12) 

Similarly, in the max-plus finite element method of Akian, 
Gaubert and Lakhoua [AGL08], the approximation is 
computed recursively by 

v\^P-^'Py,M^V) > 

where at each step, we also have a (dual) min-plus linear 
space generated by finitely many test-functions. Then, 
it follows from [AGL08] that the sup-norm projection error 
cannot be expected to be smaller than 

\\v{-,t-x)-Py^,,v{-,t-x)\\^ + \\v{-,t-x)-P^'v[-,t-T)\\^ . 

Moreover, it is shown in [AGL08] that this estimate is tight, 
meaning that the total error of the method, ||v(-,r) — Vy^||oo 
is of order at most times the previous sum, up to a term 
depending on the quality of approximation of [w] by S-^ [w] . 

Finally, in the curse of dimensionality method of McE- 
neaney [McE07], the basis functions are quadratic forms and 
the semi-group Sx is approximated by the pointwise supre- 
mum of semi-groups associated to linear-quadratic control 
problems 

5.[(/)] = v;;f=i5;"w . (13) 

The number of basis functions of the linear combination 
grows by a factor of M at each step. Then, the accuracy 
of the method is still limited by the projection error, which 
arises when pruning the less useful basis functions. 

III. Curse of dimensionality for semiconvex 

BASED approximations 

As discussed above, the main source of inaccuracy of max- 
plus basis methods is the projection error ( [T2] i. In this section, 
we give an asymptotic estimate of the optimal projection 
error as the number of basis functions tends to infinity, in 
the special case in which the space Wt consists of quadratic 
functions (|6]), as in the max-plus basis method [FMOO], or 
in the P2 type finite element method of [AGL08], [Lak07]. 

Let c e M, e > and v^(x) : M'' ^ M be a (c - e)- 
semiconvex function. It is approximated by a given number 
n of semiconvex basis functions — 5 jjcj^ + p^x: 

V/W~V/(x):== max {~^-\x\^+pJx + a{pi)) . (14) 

/= 1 n Z 

We are interested in the L\ or approximation error 

£1:= / \\if{x)--^f{x)\dx, or £00 := sup|v/(x) - ^/(jt:)! 

Jx xex 

for some suitable full dimensional compact convex subset 
X C W^. We shall compute these errors when Xj/ — v{-,t), 
and so, for reasons discussed in Section |II-C| (see (fTT)), we 
shall require that iff < \f/. We denote by d^^{\j/,c) (resp. 
5j„(V/,c)) the minimal Li (resp. L^o) approximation error 
on X of \lf{x) by ^/(x) as in ([T4]i, by xif'J the hessian matrix 
of Xj/ at point x, and by Id the identity matrix of size d. 



The next two theorems imply that whatever computation 
scheme is chosen for the coefficients fl(z,), the approximation 
error is necessarily subject to a curse of dimensionality. 

Theorem 3.1 (Li approximation error): Let c G M, e > 
and let X C M'' denote any full dimensional compact convex 
subset. If : M'' M is (c — e)-semiconvex of class 
then, there exists a constant ai > depending only on d 
such that 



Here S^{C) and 5^ (C) are respectively the minimal 
distance with respect to Hausdorff and Li metric between 
C and any circumscribed polytope with n facets, a is the 
ordinary surface area measure on dC. Moreover, we have 
the following asymptotic estimates for the constants ai and 
a2: 



ai 



d+1 



2 

nd 



{det{\l/J + cI^))d+2 dx) as n 



d±2 
d 



as c/ — > oo, 



where tJ^/ is estimated as [Rog64]: 



Theorem 3.2 (L^ approximation error): Let c, £, X and -j-^ < < dlogd + dlog(logd) +5d, with ■ 



\l/{x) be as in Theorem 3.1 Then there exists a constant 
a2> depending only on d such that 



02 

2 

nd 



(det(v/; 



-clj)) ^dx\ as n 



The proof of these theorems is sketched in the next sec- 
tion, it builds on analogous methods and results of Gru- 
ber [Gru93], [Gru07], concerning the approximation of 
smooth strictly convex bodies by circumscribed polytopes, 
the constants a\ and a2, which grow slowly with d, akeady 
appeared there. 

The following theorem is a direct corollary: 
Theorem 3.3: Assume that := v(-,r) is the value func- 
tion, and that it is and c-semiconvex. Then, for any max- 
plus basis method providing an approximation from below of 
of the value function by a supremum of n quadratic functions 
(see ([l4ji), the L} error ei and L°° error of approximation 
of the value function are both ii(^) 

as n — )• oo. 

Besides, the estimates of Theorems |3.1| and |3.2| confirm 
that when n is sufficiently large, it is more interesting to 
choose the smallest c such that ^f{x) + is convex. Note 
also that the integral term can be small if the Hessian of y/ 
is nearly constant and close to — c/^/ (attenuation of the curse 
of dimensionality). 

IV. Sketch of proof of TheoremsI3.1IandI3.2I 



In this section, we sketch the proof of Theorems 3.1 
and |3.2| Let c, e, X and satisfy the conditions of 



these theorems. Note that approximating i// as in ([14]) is 
equivalent to approximating the strongly convex function 
^f{x) + by n of its affine minorants. Hence, it suffices 
to prove Theorems 3.1 and 3.2 when c = 0, i// is strongly 
convex, and is approximated by the supremum of n of its 
affine minorants. We denote by V\// the gradient operator 
of ^, by X the point (jc, (//(x)) e M''+^ and by v4'(') the 
quadratic form determined by the hessian matrix Let 
us first recall some results on the approximation of convex 
bodies by circumscribed polytopes. 

Gruber proved in [Gru93], [Gru07] that for any convex 
body C C W^^^ with "lo^ boundary and positive Gaussian 
curvature Kc > 0, there are two constants ai and a2 depend- 
ing only on d such that as n — > oo: 



d„"{C)^a2 



Kc{x)^da{x) 



■ 1 

n'd 



{Kc{x))^da{x) 



2 

nd 



Both of the proofs partition dC into finitely many pieces 
associated to a family of supporting planes. For each support- 
ing plane, there is a corresponding strongly convex function 
whose graph is a piece of dC. Then, the volume of the 
difference between C and a circumscribed polytope can be 
estimated by computing the Li norm of the difference of this 
strongly convex functions with some of its piecewise affine 
lower bounds. For the Hausdorff metric case, some results 
regarding the optimal covering of a manifold by geodesic 
discs are used. We next apply the same techniques to our 
problem. 

First of all, we recall the definition of the Bregman 
distance. 

Definition 4.1 ([Bre67]): For any two points x and y of 
X, the Bregman distance Dy,(-;-) from x to y, associated to 
a strongly convex and differentiable function xj/, is defined 
by 

D^{x-y) = w{x)~yfiy)-Vyf(yf{x-y) . (15) 
The Bregman distance is positive definite (D^(x,3') > and 
the equality holds if and only if x = y), but it may not be 
symmetric. 

The proofs follow essentially the same steps as for convex 
bodies. We give just an outline of the proof without going 
into details. 

To prove Theorem 1 3. 1| we need an asymptotic formula for 
optimal quantization: 

Theorem 4.1 ([Gru93]): Let J C W' be Jordan measurable 
with positive volume v{J) > 0, and q a positive definite 
quadratic form on M''. Then as m — > oo; 



f d+2 N i 1 

inf I 'mm{q{s — t)\ ds ^2a\v{J) d (det^jd — y . 

SCKd.\S\=mJj t&S 

Sketch of proof of Theorem \3.1\ for c = 
Let A > 1, for each p ^ X, there is an open convex 
neighborhood U CX such that: 



1 



xi/'^ix)<xi/'!ix)<Xxi/^ix), yueu, Vxe: 



Then for every x,y £ U, the Bregman distance is 
bounded below and above as: 

^W;ix-y)<D^{x;y)<^W';Ax-y) . (16) 

We choose finitely many points pi,p2, . . . ,Pm with respec- 
tive neighborhoods t/i,...,f/„, covering the compact set X. 
One may then dissect the integral J^[nnny^sD\i/{x;y)]dx on 



smaller pieces {7, C Ui,i — I . . . ,m} with {Jj,i = 1, . . . ,m} 
Jordan measurable. Using the asymptotic formula of Theo- 
rem 4.1 and some other arithmetic inequalities as in [Gru07], 



the theorem can be deduced. 

Sketch of proof of Theorem \3.2\ for c = 

Let X be the graph of y/{x) on X. X is a li-dimensional 
(Riemannian) manifold of class with metric of class 
For each x,y G X, one may define the Riemannian metric 
between x and y, Y{x,y), by: 

7(i,y)=inf{ \ifl^^^^{u{t))^ dt\u{t) eC\u{0) =x,u{l) =y}. 



To prove Theorem |3.2| we need an asymptotic formula on 
the minimum covering. The next lemma is a special case of 
Lemma 1 in [Gru93]: 

Lemma 4.1 (Compare with Lemma 1 in [Gru93]): For 
p > 0, let n(X,p) be the minimum number of discs of 
radius p with respect to the Riemannian metric needed to 
cover X. Then: 

n{p)^ al^^l^J {det\j/!^)^dx^p'\ as p -> . (17) 

Given a similar result of minimum covering under Euclidean 
metric of Hlawka [Hla49], this lemma essentially proves the 
equivalence between the Riemannian metric on X and the 
Euclidean metric on X. Since our problem is to minimize the 
maximum radius of Bregman balls covering the manifold X, 
one last thing to be proved is the equivalence between the 
Riemannian distance and the Bregman distance. Indeed, we 
prove that: 



3M>0yx,yeX, 



Dy,{x;y) 



<M 



yx,y e Ui, Y{x,y) < disty{x,dUi] 



1 



7{x,yy<D^ix,y)<^Yix,y) 



(18) 



(19) 



Here disty{x, dUi) — mm{Y{x,y) : y E dUi}. Using the min- 
imum covering asymptotic estimation ( [T7| i, the equivalence 
between Bregman distance and Riemannian distance ([T8]l and 
(fT9]l, the desired theorem is established as in [Gru93]. 



V. Primal curse of dimensionality free method 

We consider the optimal control problem for switched lin- 
ear system studied in [McE07] (see also [MDG08], [MKIO]). 
Let ^ = {1,2,...,M}. 



y(x)=sup sup sup / L^^'i^,)- ^\(0,\^dt, (20) 



where 



L^'{x)^^x'^D^'x+{l'^'fx + a^', (21) 

^oc = {jU : [0, °o) ^ ^ : measurable}, (22) 

(23) 



w = 4°'([o,-);K'^), 

and the state dynamics are given by 



where a"' and y are such that E'" = ^a"'{a'") , Vm e 
The corresponding HJB PDE is: 

= -H = -max{//'"(x,Vy)}, (25) 



where H"' has the form 

H'"{x,p) 



\x^D"'x+\p'^irp- 



[a'^xYp + niYx + {qYp + a" 



(26) 



Under certain technical assumptions [McE07] which we will 
not repeat here, the function V is finite, it is a viscosity 
solution of ( |25] l, and it is given by V = limj-^^o St [0] where 
{St) is the Lax-Oleinik semi-group of the Hamilton- J acobi 
equation (|3j for H defined in (|25|. In [McE07], the value 
function is approximated by a max-plus sum of quadratic 
functions, and the approximated semi-group is propagated 
in a dual space. We next introduce a variant, which we 
call the primal curse of dimensionality free method: it is 
equivalent if no trimming is performed, but it avoids the use 
of dual representations. 

A. Approximate propagation 

The basis functions are allowed to be all of the quadratic 
functions smaller than the value function. Let 5"' be the 
evolution semi-group of the Hamilton-Jacobi equation Q 
for //„, defined in (|26|, m £ ^ . We approximate St\^\ by: 

5,[</)]~5,[(/)] = V„,e.#5";[0], (27) 

which for quadratic still yields a maxima of quadratics. 

B. Computation of a single semi-group operator 

The propagation of a quadratic function by S"' reduces to 
solving a differential Riccati equation (DRE). Moreover, it is 
well-known that one can recover the solution of a DRE from 
a system of Hamiltonian linear differential equations (see, 
e.g., [Rei72]). Suppose there are only quadratic terms, i.e., 
/™ = 0,/^" = 0, a'" = 0. Let (j){x) = ^x^Pqx, then S';'[(j)]{x) = 
\x^PtX, where P,^Y,X,-'^ and {X,,Y,) are the solution of: 



D" 



{A'"Y 



(28) 



X(0)=/,/, Y[Q)=Pq 

We denote by the matrix coefficient in the above linear 
system. Note that the invertibility of X{t) can be derived 
from the fact that the value function V is finite. Given a 
fixed time step T > 0, the fundamental solution exp(j2/ t) of 
the previous linear system satisfies: 



= exp(£/ t) 



(24) 



In the presence of linear or constant terms in the control 
system or in the quadratic function, the problem can be easily 
transformed into a purely quadratic one by adding a constant 
state variable. The above analysis shows that, given a fixed 
propagation time T, computing S'"[^], for every quadratic 
form 0, reduces to a matrix multiplication and an inverse 
operation, which can be done in 0{d^) incremental time. 



C. Propagation and curse of complexity 

We choose a time discretization step T, a number of steps 
K and an initial function Let ~ \/j^j^j'{x) be the 

approximation at step I. The iteration formula is given by: 

The computational growth in the space dimension is cubic, 
as shown in the above subsection. However, the number 
of quadratic forms grows by a factor of M at each itera- 
tion. This curse-of-complexity issue also occurs for the dual 
method in [McE07]. Some SDP relaxation based pruning 
method was proposed in [MDG08] to reduce the number 
of quadratic forms. We next discuss improvements of this 
pruning methods, still partly SDP based, but now exploiting 
the combinatorial nature of the problem. 

VI. Reduction of pruning to ^-center and 

^-MEDIAN PROBLEMS FOR A BREGMAN TYPE DISTANCE 

A. Formulation of the pruning problem 

We first give a general formulation for the pruning problem 
appearing in max-plus basis methods. Let F = {1,2, ... ,«/■}. 
Let 0(x):M''->Mbe a max-plus sum of «/ basis functions: 
^{x) = W jeF^i{x)- Let Q < k < nj be a fixed integer. The 
problem is to approximate (l){x) by keeping only k basis 
functions. To measure the approximation error, we introduce 
a Bregman type distance dist0(x;j) between each point x G 
and each basis function 0y(-)> such that: 

VjceM'', 3joeF, s.t. dist^(x;;o) = ; 

Vx e R'',i, j e F, dist0 (x; /) < dist^ (x; ;) 4^ ^j{x) < 0,-(x) . 

In other words, the distance dist||,(x;y) measures the loss 
at point X caused when approximating 0(-) by For 
example, the simplest choice is to let dist^ (x; j) = (/) (x) — 
(^j(x). Consider a compact set X C M'' on which we measure 
the loss. One may minimize the total loss {L\ metric) or the 
maximal loss (Loc metric) on X. 

1} L\ metric and k-median problem: 

5/(0)= min / [mindistA(x; /)lt/x . (29) 

scF,\s\=kJx jes ^ 



2) Loo metric and k-center problem: 



5r(0) = min max[mindistrf,(x; /)1 

scF.\s\=kxex ^ jes ^ 



(30) 



We recognize in ( [29) and ( |30| l the classical ^-median 
and the A:-center facility location problem with continuous 
demand area and discrete service points. The facility location 
problem, discrete or continuous, is known to be A^f-hard 
even with euclidean distance. Besides, we remark that a 
subproblem of Problem p9l ) is the volume computation for 
polytopes, which is known to be #P-hard. To the best of our 
knowledge, the only few references that discuss this general 
class of location problem replace the continuous demand 
with a discrete one with large number of points, see [DD97]. 
In the following, we consider a specific case and propose a 
method based on SDP relaxation to generate discrete points. 



B. Pruning methods 

We assume that all basis functions are quadratic: 0;(x) = 
^x^AjX + b^jX+ \cjyj G F . We normalize the distance func- 
tion as in [MDG08], i.e., dist^(x;;) = (0(x) - 0^(x))/(l + 

|X|2) . 

1) 'Sort upper bound' [MDG08]: The first method was 
introduced in [MDG08]. Roughly speaking, to each basis 
function 0j(x) we associate an importance metric : 



- maxmin((/)j(x) — (pAx))/ (1 + |x|' 



(31) 



Then Vj is the normalized error caused by pruning the 
function 0j(x). In some sense the bigger Vj is, the more 
useful the function 0y(x) is. In particular, when Vj < the 
function (j)j{x) is dominated by the others and it can be 
pruned without generating any approximation error Let 



2^ b, 



The problem ( |3T| i is equivalent to: 

Vj= max {v:3'i^O;||3'|Hl;/e;V>v,VjVj}. (32) 

vgR;>-gM''+' 

This nonconvex QCQP(quadratically constrained quadratic 
program) [BV04] has its SDP relaxation given by: 

Yn>0; Tr(F) = l; F^y/; 
Tr(ye/)>V, "ij^i. 



V j= max 

— ".V^Q 



yet 



(33) 



Then v j is an upper bound of the importance metric Vj. 
Finally the sort upper bound method consist in sorting all 
the upper bounds {v jj G F} and picking up the k first ones. 

2) 'Sort lower bound': The SDP relaxation (|33|) provides 
not only an upper bound on the importance metric but also 
a rather simple way to generate feasible solutions. 

Suppose {Y,y) is a solution of program ( (33] l. We use the 
randomization technique [FerOO] to get feasible points: we 
pick y as a Gaussian random variable with y ^ jVij,Y — 
yy^). Then over this distribution, the constraints in ( (32] l are 
satisfied on average. By sampling y a sufficient number of 
times, we get a y such that the inequality constraints in 
( |32j l are all satisfied. Then, setting x — {yi/yx, ■ ■ ■ iyd+\/yiY 
provides a lower bound of pT| ). The proposed procedure 
provides in practice a good lower bound, although there is 
no theoretical guarantee in the present generality. 

Then, the sort lower bound method proceeds as follows. 
Fixing an integer N > 0, for each basis function 0^, we 
resolve the SDP program ( |33| l, get feasible points x e Mf' 
using the above randomization technique and put them into 
a set X' . At the end we get a discrete set X' and for each 
basis function 0;(x) we calculate its lower bound by: 

xeX J f=j 

Finally we sort all of the lower bounds {y_j,j G F} and keep 
the k first ones. 

Following the above randomization technique, we get a 
discrete set X' which in some sense reflect rather well the 



importance of each basis function. We replace the compact 
set X by this discrete set X' and seek to minimize the total 
loss on X'. This gives the discrete A:-median problem: 



5~ min V [mindist^fx; / 

SCF.\S\=k .-^^ 



(34) 



This central problem in combinatorial optimization has seen 
a succession of papers designing approximations algorithms. 
Our two last pruning methods are merely two heuristics for 
the ^-median problem (|34]). 

3) 'J-V facility location': Lin and Vitter [LV92] proved 
that the constant factor approximation for general ^-median 
problem is A^P-hard. For metric distance, Jain and Vazi- 
rani [VazOl] proposed a primal-dual 6-approximation al- 
gorithm. This algorithm is interesting not only due to its 
constant factor, but also because it is combinatorial (there is 
no need to solve a linear program). 

4) 'greedy facility location': The fourth method is the 
greedy heuristic. Remember that the function to be mini- 
mized in the facility location problem is supermodular, which 
implies that the greedy heuristic has a bound estimate (even 
without the triangular inequality on the distance function). 
Let 5g be the value of a particular solution constructed by the 
greedy heuristic, then we have [NWF78]: 5^ < (1 — a'^)5 
a*' {max j^F Hxex' '^^^^ {^'^ J)) , where a 
time of the greedy heuristic is 0{km). 



k-\ 



The execution 



VII. Experimental results 
A. Problem instance 

To compare with the dual max-plus basis method, we use 
the instance of [MDG08] originating from H-infinity control, 
in which the parameters where chosen to exhibit a complex 
behavior. The state dimension and the switch number are 
both 6. The overpruning threshold is also the same: we keep 
k{i) =20 + 6/ basis functions at step /. 

Without the exact value function, we do not have a direct 
error estimation. Recall that the value function V is the 
unique viscosity solution of the following HJB equation: 







H 



-max{i/'"(jc,Vy)}, 



(35) 



where H'^ is defined in ( [26| l. The value of Hamiltoruan is 
then used to measure the approximation. 

B. Numerical results 

All of our result^ are shown along the x\-X2 axes with 
the 4 other coordinates of x set to 0. 

For comparison with [MDG08], we first take the same 
time-discretization step-size T — 0.2, the same iteration steps 
25 and the same pruning method {sort upper bound). Fig- 
ure [T] shows the value of Hamiltonian H at the end of 25 
iterations. Comparing with the error plot shown in [MDG08], 
which is in the same scale but has a peak of error of order 

'The code was mostly written in Matlab (version 7.1L0.584), calling 
YALMIP (version 3) and SeDuMi (version 1.3) for the resolution of SDP 
programs. The computation of the distance function and Jain & Vazirani's 
primal dual algoiithm were wiitten in C++. The results were obtained on a 
single core of an Intel quad core running at 2.66GHz, with 8Gb of memory. 



1 (versus 0.3 here), we see that the primal max-plus basis 
method yields a small improvement. 




Fig. 1: Backsubstitution error on the xi,X2 plane, T = 0.2, 
with sort upper bound pruning method 

Now we take smaller time-discretization step-sizes. Fig- 
ure [2] compares the four pruning methods with T = 0.1 and 
T — 0.05. They both show that the sort lower bound and the 
greedy facility location pruning method are better than the 
two others. 
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Fig. 2: Discrete L\ norm of H, T = 0.1 (top), T ~ 0.05 
(bottom) 



C. Discussion 

Our experimental results confirm that the total approxi- 
mation error comes both from the approximation error of 
the Lax-Oleinik semi-group and from the pruning error. The 
error of approximation of the semi-group can be improved 
by decreasing the discretization-time step-size T, while the 



TABLE I: CPU time 



T=0.2, K=25 


Total time 


Propagation 


SDP 


Pruning 


sort lower 


1.04h 


1.85% 


98.15% 


0.00% 


sort upper 


1.34h 


1.52% 


98.43% 


0.05% 


J-V p-d 


1.38h 


1.45% 


89.47% 


9.08% 


greedy 


1.43h 


1.63% 


97.84% 


0.53% 



pruning error depends on the pruning method. When T 
becomes small, the pruning appears to be the bottleneck. 
We introduced here new pruning methods, combining facil- 
ity location algorithms and semidefinite relaxations, which 
improve the final precision (see Figure [2] ). However, these 
pruning methods remain time-consuming (see Table new 
ideas are needed to develop more efficient methods. Our 
experiments also show that the error is of order 0{t), which 
is smaller than the bound of 0{\/t) established in [MKIO]. 
This remains to be studied theoretically. 





-2 -2 



Fig. 3: Backsubstitution error (top) and Optimal policy 
(bottom) on the xi,X2 plane, T = 0.1, with the greedy facility 
location pruning method 
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